Force mobilization and generalized isostaticity in jammed packings of frictional grains 
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We show that in slowly generated 2d packings of frictional spheres, a significant fraction of the 
friction forces lies at the Coulomb threshold — for small pressure p and friction coefficient p,, about 
half of the contacts. Interpreting these contacts as constrained leads to a generalized concept of 
isostaticity, which relates the maximal fraction of fully mobilized contacts and contact number. For 
p — > 0, our frictional packings approximately satisfy this relation over the full range of p,. This 
is in agreement with a previous conjecture that gently built packings should be marginal solids at 
jamming. In addition, the contact numbers and packing densities scale with both p and fx. 

PACS numbers: 45.70.-n, 46.65.+g, 83.80.Fg 



Models of frictionless polydisperse particles with finite- 
range repulsive forces exhibit a well-defined "jamming 
point" J in the limit that the confining pressure p goes 
to zero 0,0- In the vicinity of J on the jammed side, i.e. 
for p > 0, the average contact number, packing density, 
elastic constants, vibrational modes and response func- 
tions all show scaling behavior as a function of pressure 
|H 13- • This scaling is intimately connected to the fact 
that when point J is approached by preparing packings 
at lower and lower pressures, such packings become iso- 
static: a simple constraint counting argument for hard 
spheres in d dimensions yields that for p — > 0, the aver- 
age number of contacts per interactingparticle z equals 
the fictionless isostatic value z® so = 2d [a, Q . 

The picture which is emerging for frictional packings 
is much more diffuse, since there are now two control pa- 
rameters (p and /j,), and more importantly, packing den- 
sities and contacts numbers depend on the preparation 
method and history. This is because the Coulomb condi- 
tion for the frictional force is an inequality: it specifies, 
for each static contact, that the tangential force / t be less 
than or equal to the friction coefficient \x times the nor- 
mal force / n : |/t| < pf n . If in view of this inequality we 
treat these tangential forces as independent new degrees 
of freedom in the constraint counting, the isostatic value 
jumps from zP = 2d to z- 1 = d+ 1, and in d dimension 
frictional packings for p — > can in principle occur for 
any z in the range z? SQ = d + 1 < z < z° s [Tfl • 

In practice, however, for a given experimental || or 
numerical |9l llCllll| protocol some reproducible value z is 
found. The sudden jump of the isostatic contact number 
with p is not reflected in a jump of zj(p) = z(p,p — > 0): 
numerically, is found to vary smoothly from z? 

at small /z to some limiting value at large \x The 
large \x limit may or may not coincide with z? , and z 
is generally smaller and closer to the isostatic value the 
slower the packing is prepared |Tlj . 

As stressed by Silbert et al. and Bouchaud [l^ . 
there is a natural way in which the discontinuity in 
the isostatic contact numbers is not reflected in Zj(/j,), 
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FIG. 1: Relation between the number of fully mobilized forces 
per particle n m and contact number z. The full line indicates 
the maximum of n m and such packings are marginal, while 
below this line one finds hyperstatic stable packings. The 
data points refer to numerically obtained values of n m in 2 
dimensions, for p ~ 2 x I0 -4 (+), p ~ 5 x I0 -5 (□), p ~ 
2 x 10 -5 (x) and p ~ 5 x I0 -6 (o). n m and z behave smoothly 
as function of p 1 ^ 3 , and by extrapolation we obtain our p = 
estimate indicated by the black squares (see inset and main 
text). 



which hinges on the notion of maximizing the number 
of fully mobilized or "plastic" contacts, i.e., those at 
the Coulomb failure threshold for which m = 1, where 
to = |/t|/(/// n ) UU E3 Since at fully mobilized con- 
tacts tangential and normal forces are related, this leads 
to additional constraints in the counting arguments: In- 
troducing the number of fully mobilized contacts 
per particle in a packing with N{ interacting particles, the 
zdNi/2 force degrees of freedom should be larger than the 
total number of constraints provided by the N[d(d+ l)/2 
force and torque balance equations and the n m Ni mo- 
bilization constraints. This gives: 



<z-zr 
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From this point of view, packings with n m = z 
in fact isostatic or marginal, while packings with n m < 



Z L are 
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z — z? Q are hyperstatic (see Fig. 0) . 

In this paper, we will show that gently prepared pack- 
ings support this scenario over a surprisingly wide range 
of friction coefficients. The distribution function P(m) 
of such packings indeed naturally splits up in a peak at 
m = 1 and a broad flat part for m < 1 (Fig. 2), and these 
packings actually tend to be marginal at jamming, i.e., 
to lie close to this generalized isostaticity line in Fig. ^ 
The picture that emerges is that if we prepare the pack- 
ings sufficiently slowly, they get stuck in a marginal state. 
Such a marginal scenario also occurs in, e.g., spinglasses 
1^1 , charge density waves 0] and phase organization 

The fact that our well-equilibrated packings approach 
a well-defined limit opens up the possibility to study the 
asymptotic scaling behavior as a function of pressure and 
friction coefficient \i. We have therefore also investigated 
the effect of the applying pressure on repeatedly and gen- 
tly created packings over a whole range of friction coef- 
ficients, and find that contact numbers z and packing 
densities <f> of the packings do exhibit scaling with p and 
li. The scaling of cf> and z with p are related to the form of 
the interparticle potential and consistent with previous 
findings for the frictionless case. The scaling of z and cf> 
with [i appear to be independent of the force law — we 
have at present no good physical understanding of this 
scaling. 

Model and simulation method — We numerically build 
2d packings of iV p = 1000 polydisperse spheres that in- 
teract through 3d Hertz-Mindlin forces or through one- 
sided-linear-springs-plus-friction 0] m a square box with 
periodic boundary conditions. The data reported below 
are all for the 3c? Hertz-Mindlin forces. Following 0] 
our units are such that the mass density, the average 
particle diameter and the Young's modulus of the grains 
are 1. The Poisson ratio of the grains is taken to be 
zero, and there is no gravity. As in the packings are 
constructed by cooling an initial low density state where 
the particles have a small velocity, while slowly inflating 
the particle radii by multiplying them with a common 
scale factor r s . This factor is determined by solving the 
damped equation r" = —Auj^r' s — [p(*> r s)/P — l] r s) 
where cj ~ 6 * 10 -2 , p(t,r s ) is the instant value of the 
pressure and p the target pressure. This ensures a very 
gentle equilibration of the packings. In our analysis of 
forces and contact numbers, we always take out rattlers, 
by considering contact forces less than 10~ 3 times the av- 
erage force broken and removing particles with less than 
two contacts. For each packing, we determine the total 
number of contacts N c and the total number of inter- 
acting particles iVj (the total number of particles minus 
the rattlers) — z = 2N C /Ni. For each value of p and 
/i € [10 _3 ,10 3 ], 30 realizations have been constructed 
with a polydispersity of 20%. We occasionally checked 
that taking 60 realizations, a different polydispersity or 
different damping parameters lead to similar results. In 




FIG. 2: Mobilization at p = 2 x 10" 5 . (a) Scatter plots 
of ft/p versus / n for three packings at p = 0.001,0.32 and 
1. The probability density of normalized tangential (ft/p) 
and normal (/„) forces exhibits a singularity on the Coulomb 
cone for small p, which rapidly diminishes for larger p (all 
forces are normalized so that {/„) = 1). (b) The cumula- 
tive distribution of the mobilization C(m) = f" 1 dm' P(m') 
exhibits a clear jump near m = 1 — data shown here is for 

fi = io°,io-°' 5 ,io- 1 ,...,io- 3 . 

comparison with other simulations where the particles 
settled under gravity or were quenched rapidly 
our algorithm prepares the packings more gently, in the 
sense that it results in low packing densities and coordi- 
nation numbers. 

The density n m of fully mobilized contacts — The joint 
probability distribution of the normal and frictional con- 
tact forces clearly show that for small fj,, a substantial 
amount of forces lie on the Coulomb cone, i.e., have 
to — 1, while for larger /i the fraction of fully mobilized 
contacts diminishes (Fig. 2a). A priori it would appear to 
be difficult to determine numerically whether a contact 
is fully mobilized with to = 1 or elastic (non-mobilized) 
with to < 1, but as Fig. 2b shows, the cumulative dis- 
tribution C(m) = J dm'P(m') exhibits a clear jump at 
to = 1. The value of n m equals z/2 [1 — C(m— ► 1)], and 
we find that for small friction about half of the contacts 
(one contact per particle) is at the Coulomb treshold! 
Especially for small fi, C(m) is linear in to, which means 
that the distribution of non-mobilized forces is flat — in 
other words, non-mobilized contacts are not biased to- 
wards higher contact numbers. 

Our estimates for n m and z for p — > and a range of /x 
lie very close to the generalized isostaticity line (Fig.^l. 
Note that we have extrapolated contact numbers and n m 
to estimate the zero pressure limit (see the inset of Fig. 1 
and Fig. 3). The close proximity of n m and z to the 
marginal line presents, to our knowledge, the strongest 
support to date for the marginal solid scenario described 
above: when frictional packings are sufficiently gently 
prepared, they form a marginally stable jammed solid 
which in a generalized sense is an isostatic solid. We 
expect that the deviations from the generalized isostatic- 
ity will be larger the faster the granular particles are 
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FIG. 3: Variation of contact numbers z and packing density <j> 
as function of pressure p and friction coefficient \x. Errorbars 
are smaller then the symbol size, (a-c) The variation of the 
contact number z, the packing density including rattlers 4>+r, 
and the packing density excluding rattlers (fi-n as a function 
of fi. Symbols indicate data at pressures p~4x 10 _3 ( V), 5 x 
l(T 4 (o),2 x I0~ 4 (+),5 x I0- 5 (D),2 x I(T 5 (x),5 x f(T 6 (o). 
Based on the extrapolation illustrated in panels (d-f), we also 
show the estimated values at p — (■). Even though <f>+n 
and 0_r differ substantially, their variation with fi appears 
very similar, (d-f) z scales as p 1 ^ 3 and as p 2 ^ 3 , which 
allows us to extrapolate to zero pressure. Surprisingly, the 
packing density 0_r does not scale convincingly with p 2 ^ 3 , 
but rather as p 1 ^ 3 . Symbols are as in panel a-c. 



compressed or quenched; earlier simulations already give 
indications for this |l(J, [llj . 

Scaling behavior of z and <f> — Since our packings for 
small p approach the generalized isostaticity line, one 
may wonder how contact number and packing density 
<f> change when moving away or along this line. Since the 
number of rattlers is strongly dependent on the pressure 
p and on the friction coefficient /i, we have found it il- 
luminating to study both the density with the rattlers 
excluded and included, <p-n and </>+r, respectively. Note 
that for small pressure and small friction about 4% of the 
particles are rattlers, which rises to 12% for large values 
of the friction. The results of our analysis are shown in 
Fig. 3a-c. As a function of /Lt, the overall variation of z in 
Fig. 3a is very similar to results obtained by contact dy- 
namics 0, and again the density variations in Fig. 3bc 
mimic that of z. As a function of p, our data is con- 
sistent with the scaling relation z(fi,p) — z(fi, 0) ~ p 1 / 3 
(Fig. 3d). This allows us to extrapolate with confidence 
to zero pressures, giving z(fi < 1,0) = 3.98 ± 0.02 and 
z(fi ^> 1,0) = 3.00 ± 0.02, which are close to the fric- 
tionless and frictional isostatic bounds, Z;° so = 4 and 
z \so = 3, respectively. For the whole range of [i we find 
that the change in density including rattlers scales as 
<f> + R(fj,,p) - 0+r(a*, 0) ~ p 2/3 (Fig. Et). This is consis- 
tent with the scaling of the density in frictionless pack- 



ings upon compressing a given packing [2| , and with the 
variation K ~ (d(/> + R/dp) _1 ~ p 1 / 3 of the compression 
modulus K with pressure 0, El ■ Interestingly, the den- 
sity excluding rattlers, 0_r appears to vary instead as 
p 1 / 3 (Fig. 3f). 

For our Hertz-Mindling forces, the p 1 / 3 scaling for z is 
consistent with the scaling z — z® so ~ y/S observed also 
for frictionless particles 0, 0], where S is the typical 
dimensionless overlap of the particles. We have checked 
that our results do only trivially depend on the details of 
the force law: for one-sided harmonic springs the z and 
4> scale as function of p 1 / 2 (not shown). The fact that 
z scales with p similarly as for frictionless systems was 
seen in some studies [llj but not in others 40j. Both the 
presence of this scaling and fact that our packings reach 
the generalized isostaticity line for p — > may be related 
to our very slow rate of equilibration. 

From the zero pressure extrapolations discussed above, 
we can study the variation of the contact number and 
densities at jamming. The results of this analysis are 
summarized in Fig. 4, with details given in the figure cap- 
tion. In particular we find z(/z,0) to decrease for small 
(i as /i°- 7±01 . That indeed z decreases rapidly with p is 
also clear from the 3d data of 01 > which appear to fit a 
powerlaw behavior Az ~ p 05 reasonably well. Whether 
the density changes for small /i with a nontrivial expo- 
nent different from 1 is less clear from our data. We can 
not draw any firm conclusion from our data regarding the 
functional p-dependence for large friction but the varia- 
tion of contact number with density appears is consistent 
with an exponent of 1.7. Similar scalings are obtained for 
linear instead of Hertzian contact laws. 

Summary and Outlook — Our results substantiate the 
scenario that when a packing is gently prepared, it gets 
jammed in a (near) marginal state, where enough con- 
tacts get stuck at the Coulomb failure threshold to make 
the packing a marginal solid. Note that this is different 
from what engineers refer to as "incipient failure every- 
where" — the idea that one can deal with the Coulomb 
inequality by turning it into an equality for all contacts 
[l8| . Our results here show that this overestimates the 
number of fully mobilized contacts. Our results suggest 
a lower boundary for the contact number, and possibly 
for the packing densities too, that can be obtained for 
finite /z, whereas naive counting would suggest that d- 
dimensional packings with any contact number between 
d + 1 and 2c? could arise. 

An immediate implication of our results is that the re- 
sponse properties of such gently prepared packings will 
have a strong tendency to show nonlinear response, de- 
pending very sensitively on the behavior of the plastic 
contacts: if these remain fixed at the Coulomb threshold, 
the fact that these packings are near isostaticity will give 
many low-frequency modes and will make these packings 
very soft. If these contacts yield, however, irreversibility 
effects will dominate. 
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FIG. 4: Scaling of the zero pressure, extrapolated, contact 
numbers and packing densities with the friction coefficient 
fi. The extrapolated values at zero (infinite) friction are la- 
belled as 0,0 (0, oo). (a-c) When /j, — > and p — > 0, z 
approaches ~ 3.975[20], while </>+r approaches </>+ti ~ 
0.8395[5]. For finite but small n, z and 4>+r appear to scale 
as: (a) (z°'° - z) ~ n ^ and (b) (<^° - +R ) ~ M - 77 [ 10 l . 
(c) The contact number and packing deviate similarly from 
this scaling when n approaches one, and so (z°'° — z) ~ 
(<£°'° - 0+r) - 911101 . (d) In the limit of infinite friction and 
zero pressure, z approaches z 00,0 = 3.00[2], while 4>+r ap- 
proaches = 0.758 [10]. The deviations from this limiting 
values also appear to be related by a scaling relation of the 
form ( Z -^V(^-C) L7f21 - 

The contact numbers and densities that characterize 
gently prepared packings show various nontrivial scaling 
relations as a function of jjL and p. The scaling of z ~ p 1 ^ 3 
and 4>+r ~ p 2 ^ 3 with p are similar to those found for 
frictionless Hertzian packings - but these scalings seem 
to work equally well over the whole range of //. The scal- 
ing of 4>-r is more puzzling. It is very well possible that 
the asymptotic behavior for very small p crosses over to 
the familiar p 2 ! 3 behavior, but we can not access this 
regime at present. In addition, for 3d packings the frac- 
tion of rattlers may be smaller than for 2d, so that there 
we expect less of this effect. Nevertheless, the question 
whether one should include or exclude rattlers is subtle 
— see also 

The scaling of z and <P+r with fi is new and presently 
not understood, but may give indirect evidence for strong 
correlations between the tangential forces. Suppose we 
think of the tangential forces ft as small randomly point- 
ing perturbations of the net forces on the particles for 
/i <C 1. In a domain of linear scale L, these tangential 
forces add up to a total force of order [if n L d / 2 . This 
is comparable to the normal force scale / n on a scale 
L M ~ [i~ 2 l d . It might therefore be natural to suppose 
that on this scale the tangential forces allow to reduce z 
by replacing a single contact. Since AzlA — 0(1), this 
would suggest Az ~ /i 2 , in strong contrast to the data. 
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